! Viscous Fluxes Calculation
subroutine CalcViscousFluxes()
use Data_Mesh, only: nCells,nFaces,nPoints, nInternalFaces, &
                     CellViscousStencil,FaceViscousStencil,PointViscousStencil, &
                     meshFaceOwner, meshFaceNeighbour, &
                     meshSf, meshFaces, meshVPoints, &
                     nBoundaryFields, boundaryMeshStart, boundaryMeshSize, &
                     opposingOwnerFaceLabel, meshCells, meshV, meshFaceCenters
use DataCells, only: CurrentVelocities
use DataFaces, only: FaceVarsNew, ViscousFluxes
use Data_Cfoam, only: GradU, btype, sumViscousCellFluxes
use DataPerfectGas, only: TemperaturePRho, DynamicViscosity, TotalBulkViscosity
use DataTime, only: TimeCell, dtCell, TimeEpsilon, dtFace, dt2Cell, CurrentTimeCell, &
                    NextTimeCell, NextTimeLocal, nNextTimeCell,nNextTimeLocal, &
                    nNextTimeCell1
implicit none

    integer :: i,j,k, bmeshs,bmeshf, m, iBound
    integer :: iPoint, iFace, iCell, CellOwner,CellNeighbour, iSide, numFace, iCellOwner
    real(8) :: Rho,Temp, DynVisc, BulkViscTotal
    real(8) :: deltaU, faceGradUij, faceGradU(3,3)
    real(8) :: stress(3,3), divVel, ViscFlux(3)
    real(8) :: minVal1,maxVal1, minVal2,maxVal2, minVal3,maxVal3
    real(8) :: grad(3,3)
    
    real(8),allocatable :: faceGradU_view(:,:,:) ! sumViscousCellFluxes(:,:),

    call DefineCurrentGlobalViscousStencil()
    call CalcExtrapolatedVelocities()
    
    allocate(faceGradU_view(3,3,nFaces))

    do iPoint = 1,nPoints
!        if(PointViscousStencil(iPoint))then
            do j = 1,3
                do i = 1,3
                    gradU(i,j,iPoint) = 0d0
                end do
            end do
!        end if
    end do

    ! compute GradU in internal points (nodes)
    do iFace = 1,nInternalFaces ! nFaces
        if(FaceViscousStencil(iFace).eqv..false.) cycle

    	do j=1,3
    	    CellOwner = meshFaceOwner(iFace)
    	    if(iFace.gt.nInternalFaces)then
    	        deltaU = CurrentVelocities(j,CellOwner) - FaceVarsNew(j+1,1,iFace)
    	    else
    	        CellNeighbour = meshFaceNeighbour(iFace)
			    deltaU = CurrentVelocities(j,CellOwner)-CurrentVelocities(j,CellNeighbour)
			end if
			do i=1,3
				faceGradUij = meshSf(i,iFace) * deltaU
				do k = 1,4
					gradU(i,j,meshFaces(k,iFace)) = gradU(i,j,meshFaces(k,iFace)) + faceGradUij ! m2 * m/s
				end do
			end do
		end do
    end do
    
    ! compute GradU in boundary points (nodes)
	do iBound=1,nBoundaryFields
        bmeshs = 1 + boundaryMeshStart(iBound)
        bmeshf = boundaryMeshStart(iBound) + boundaryMeshSize(iBound)

		if(btype(iBound).eq.3)then
		    
			do iFace = bmeshs,bmeshf
			    CellOwner = meshFaceOwner(iFace)
			    m = opposingOwnerFaceLabel(iFace)
                if(meshFaceOwner(m).eq.CellOwner)then
                    iSide = 1
                else
                    iSide = 2
                end if
				do j = 1,3
					deltaU = (4d0*CurrentVelocities(j,CellOwner) - FaceVarsNew(j+1,iSide,m)) * 0.5d0 ! / 3d0

			        do i=1,3
				        faceGradUij = meshSf(i,iFace) * deltaU
				        do k = 1,4
					        gradU(i,j,meshFaces(k,iFace)) = gradU(i,j,meshFaces(k,iFace)) + faceGradUij ! m2 * m/s
				        end do
			        end do
				end do
			end do
		end if

		if(btype(iBound).eq.4)then
            cycle
		end if

		if(btype(iBound).eq.1.or.btype(iBound).eq.2)then ! INLET, OUTLET
			do iFace = bmeshs,bmeshf
			    CellOwner = meshFaceOwner(iFace)
			    m = opposingOwnerFaceLabel(iFace)
                if(meshFaceOwner(m).eq.CellOwner)then
                    iSide = 1
                else
                    iSide = 2
                end if
				do j = 1,3
					deltaU = CurrentVelocities(j,CellOwner) - FaceVarsNew(j+1,1,iFace)
!					deltaU = (4d0*CurrentVelocities(j,CellOwner) - FaceVarsNew(j+1,iSide,m)) * 0.5d0

			        do i=1,3
				        faceGradUij = meshSf(i,iFace) * deltaU
				        do k = 1,4
					        gradU(i,j,meshFaces(k,iFace)) = gradU(i,j,meshFaces(k,iFace)) + faceGradUij ! m2 * m/s
				        end do
			        end do
				end do
			end do
		end if
	end do
    ! ************************************************************************************************* !

	do iPoint = 1,nPoints
	    if(PointViscousStencil(iPoint).eqv..false.) cycle

		do i=1,3
			do j=1,3
				grad(i,j)=gradU(i,j,iPoint)/meshVPoints(iPoint) !  ! 1 / s
			end do
		end do

		do j=1,3
			do i=1,3
				gradU(i,j,iPoint) = grad(i,j) + grad(j,i) !  ! 1 / s
			end do
		end do
    end do

    do iFace = 1,nInternalFaces
        CellOwner = meshFaceOwner(iFace)
        CellNeighbour = meshFaceNeighbour(iFace)
!        if(nNextTimeCell1(CellOwner).eqv..false..or.nNextTimeCell1(CellNeighbour).eqv..false.)then
        if(nNextTimeCell1(CellOwner).or.nNextTimeCell1(CellNeighbour)) cycle

            faceGradU = 0d0

            do k = 1,4
                do j = 1,3
                    do i = 1,3
                        faceGradU(i,j) = faceGradU(i,j) + gradU(i,j,meshFaces(k,iFace))
                    end do
                end do
            end do
            
            do j = 1,3
                do i = 1,3
                    faceGradU(i,j) = 0.25d0 * faceGradU(i,j)
                end do
            end do
    
            faceGradU_view(1:3,1:3,iFace) = faceGradU(1:3,1:3)
            faceGradU_view(1,1,iFace) = 0.5/1.5e-5 * (12*1.5e-5/5d-3**2*0.01)*(0.005-2*meshFaceCenters(3,iFace))

            if(nNextTimeCell1(CellOwner).eqv..false.)then
                Rho = FaceVarsNew(1,1,iFace)
                Temp = TemperaturePRho(FaceVarsNew(5,1,iFace),Rho)
            else
                Rho = FaceVarsNew(1,2,iFace)
                Temp = TemperaturePRho(FaceVarsNew(5,2,iFace),Rho)
            end if

            divVel = faceGradU(1,1) + faceGradU(2,2) + faceGradU(3,3)

            DynVisc = DynamicViscosity(Temp,Rho)
		    BulkViscTotal = TotalBulkViscosity(DynVisc)

            stress = 0d0
            
            do j = 1,3
                do i = 1,3
                    stress(i,j) = stress(i,j) + DynVisc * (faceGradU(i,j)) ! + faceGradU(j,i))
                    if(i.eq.j)then
                        stress(i,j) = stress(i,j) + 2d0 * BulkViscTotal * divVel
                    end if
                end do
            end do

            do i = 1,3
                ViscFlux(i) = 0d0
                do j = 1,3
                    ViscFlux(i) = ViscFlux(i) + stress(i,j) * meshSf(j,iFace)
                end do
            end do

!meshSf(1,k)*stress(1,1:3)+meshSf(2,k)*stress(2,1:3)+meshSf(3,k)*stress(3,1:3)

            faceGradU_view(1,1,iFace) = faceGradU_view(1,1,iFace) * 1.5e-5 * sqrt(sum(meshSf(1:3,iFace)*meshSf(1:3,iFace)))

            if(nNextTimeCell1(CellOwner).eqv..false.)then
                do i = 1,3
                    ViscousFluxes(i,1,iFace) = -ViscFlux(i)
                end do
            end if

            if(nNextTimeCell1(CellNeighbour).eqv..false.)then
                do i = 1,3
                    ViscousFluxes(i,2,iFace) = ViscFlux(i)
                end do
            end if

!        end if
    end do
    
    ! external faces
    
    do iFace = nInternalFaces+1,nFaces
        CellOwner = meshFaceOwner(iFace)
!        if(nNextTimeCell1(CellOwner).eqv..false.)then
        if(nNextTimeCell1(CellOwner)) cycle

            do j = 1,3
                do i = 1,3
                    faceGradU(i,j) = 0d0
                    do k = 1,4
                        faceGradU(i,j) = faceGradU(i,j) + gradU(i,j,meshFaces(k,iFace))
                    end do
                    faceGradU(i,j) = 0.25d0 * faceGradU(i,j)
                end do
            end do
            faceGradU_view(1:3,1:3,iFace) = faceGradU(1:3,1:3)
            
            faceGradU_view(1,1,iFace) = 0.5/1.5e-5 * (12*1.5e-5/5d-3**2*0.01)*(0.005-2*meshFaceCenters(3,iFace))
            
            Rho = FaceVarsNew(1,1,iFace)
            Temp = TemperaturePRho(FaceVarsNew(5,1,iFace),Rho)

            divVel = faceGradU(1,1) + faceGradU(2,2) + faceGradU(3,3)

            DynVisc = DynamicViscosity(Temp,Rho)
		    BulkViscTotal = TotalBulkViscosity(DynVisc)

            stress = 0d0
            
            do j = 1,3
                do i = 1,3
                    stress(i,j) = stress(i,j) + DynVisc * faceGradU(i,j) !+ faceGradU(j,i))
                    if(i.eq.j)then
                        stress(i,j) = stress(i,j) + 2d0 * BulkViscTotal * divVel
                    end if
                end do
            end do

            do i = 1,3
                ViscFlux(i) = 0d0
                do j = 1,3
                    ViscFlux(i) = ViscFlux(i) + stress(i,j) * meshSf(j,iFace)
                end do
            end do

            do i = 1,3
                ViscousFluxes(i,1,iFace) = -ViscFlux(i)
            end do
!        end if
    end do
    
    minVal1 = minVal(ViscousFluxes(1,1,1:nFaces))
    maxVal1 = maxVal(ViscousFluxes(1,1,1:nFaces))

    minVal2 = minVal(ViscousFluxes(2,1,1:nFaces))
    maxVal2 = maxVal(ViscousFluxes(2,1,1:nFaces))
    
    minVal3 = minVal(ViscousFluxes(3,1,1:nFaces))
    maxVal3 = maxVal(ViscousFluxes(3,1,1:nFaces))

!    allocate(sumViscousCellFluxes(3,nCells))
    sumViscousCellFluxes = 0d0
    
    do i = 1,nCells
        do numFace = 1,6
            iFace = meshCells(numFace,i)
            iCellOwner = meshFaceOwner(iFace)
            if(iCellOwner.eq.i)then
                iSide = 1
            else
                iSide = 2
            end if
            sumViscousCellFluxes(1:3,i) = sumViscousCellFluxes(1:3,i) + ViscousFluxes(1:3,iSide,iFace)
        end do
!        sumViscousCellFluxes(1:3,i) = sumViscousCellFluxes(1:3,i) / meshV(i)
    end do

end subroutine CalcViscousFluxes


subroutine DefineCurrentGlobalViscousStencil()
use Data_Mesh, only: nCells, nFaces, nPoints, &
                     nCell2NeigCells, Cell2NeigCells, &
                     nFace2NeigCells, Face2NeigCells, &
                     nPoint2NeigCells, Point2NeigCells, &
                     CellViscousStencil, FaceViscousStencil, PointViscousStencil
use DataTime, only: nNextTimeCell1
implicit none

    integer :: iCell,jCell,iFace,iPoint
    integer :: i

    ! compute cell global viscous stencil
    do iCell = 1,nCells
        CellViscousStencil(iCell) = .false.
    end do

    do iCell = 1,nCells
        do i = 1,nCell2NeigCells(iCell)
            jCell = Cell2NeigCells(i,iCell)
            if(nNextTimeCell1(jCell)) cycle
            CellViscousStencil(iCell) = .true.
            exit
        end do
    end do

    ! compute face global viscous stencil
    do iFace = 1,nFaces
        FaceViscousStencil(iFace) = .false.
    end do

    do iFace = 1,nFaces
        do i = 1,nFace2NeigCells(iFace)
            jCell = Face2NeigCells(i,iFace)
            if(nNextTimeCell1(jCell))cycle
            FaceViscousStencil(iFace) = .true.
            exit
        end do
    end do
    
    
    ! compute point global viscous stencil
    do iPoint = 1,nPoints
        PointViscousStencil(iPoint) = .false.
    end do

    do iPoint = 1,nPoints
        do i = 1,nPoint2NeigCells(iPoint)
            jCell = Point2NeigCells(i,iPoint)
            if(nNextTimeCell1(jCell))cycle
            PointViscousStencil(iPoint) = .true.
            exit
        end do
    end do

end subroutine DefineCurrentGlobalViscousStencil

subroutine CalcExtrapolatedVelocities()
use Data_mesh, only: nCells, CellViscousStencil
use DataTime, only: nNextTimeCell, nNextTimeLocal, dnTimeCell
use DataCells, only: CurrentVelocities
use Data_Cfoam, only: Unew, U
implicit none

    integer :: iCell, i
    real(8) :: alfa

    do iCell = 1,nCells
        if(CellViscousStencil(iCell).eqv..false.) cycle

        do i = 1,3
            alfa = 2.0d0 * (nNextTimeLocal + dnTimeCell(iCell) - nNextTimeCell(iCell)) * (real(dnTimeCell(iCell)))**(-1)
! linear extrapolation
!            alfa = (nNextTimeLocal + dnTimeCell(iCell) - nNextTimeCell(iCell)) * dnTimeCell(iCell)**(-1)
!            CurrentVelocities(i,iCell) = U(i,iCell) + (Unew(i,iCell)-U(i,iCell)) * alfa

! Adam's method
! U_n+1 = U_n + tau_n+1 * F_n = U_n + tau_n+1 / tau_n+1 * (U_n - U_n-1)
!            alfa = 0d0 ! (nNextTimeLocal - nNextTimeCell(iCell)) * dnTimeCell(iCell)**(-1)
            CurrentVelocities(i,iCell) = U(i,iCell) + (Unew(i,iCell)-U(i,iCell)) * alfa
!            CurrentVelocities(i,iCell) = Unew(i,iCell) + (Unew(i,iCell)-U(i,iCell)) * alfa
!            CurrentVelocities(i,iCell) = Unew(i,iCell)
        end do
    end do

end subroutine CalcExtrapolatedVelocities


!! Viscous Fluxes Calculation
!subroutine CalcViscousFluxes()
!use Data_Mesh, only: nCells,nFaces,nPoints, nInternalFaces, &
!                     CellViscousStencil,FaceViscousStencil,PointViscousStencil, &
!                     meshFaceOwner, meshFaceNeighbour, &
!                     meshSf, meshFaces, meshVPoints
!use DataCells, only: CurrentVelocities
!use DataFaces, only: FaceVarsNew, ViscousFluxes
!use Data_Cfoam, only: GradU
!use DataPerfectGas, only: TemperaturePRho, DynamicViscosity, TotalBulkViscosity
!use DataTime, only: TimeCell, dtCell, TimeEpsilon, dtFace, dt2Cell, CurrentTimeCell, &
!                    NextTimeCell, NextTimeLocal, nNextTimeCell,nNextTimeLocal, &
!                    nNextTimeCell1
!implicit none
!
!    integer :: i,j,k
!    integer :: iPoint, iFace, iCell, CellOwner,CellNeighbour
!    real(8) :: Rho,Temp, DynVisc, BulkViscTotal
!    real(8) :: deltaU, faceGradUij, faceGradU(3,3)
!    real(8) :: stress(3,3), divVel, ViscFlux(3)
!
!    call DefineCurrentGlobalViscousStencil()
!    call CalcExtrapolatedVelocities()
!
!    do iPoint = 1,nPoints
!        if(PointViscousStencil(iPoint))then
!            do j = 1,3
!                do i = 1,3
!                    gradU(i,j,iPoint) = 0d0
!                end do
!            end do
!        end if
!    end do
!
!    ! compute GradU in points (nodes)
!    do iFace = 1,nFaces
!        if(FaceViscousStencil(iFace).eqv..false.) cycle
!
!    	do j=1,3
!    	    CellOwner = meshFaceOwner(iFace)
!    	    if(iFace.gt.nInternalFaces)then
!    	        deltaU = CurrentVelocities(j,CellOwner) - FaceVarsNew(j+1,1,iFace)
!    	    else
!    	        CellNeighbour = meshFaceNeighbour(iFace)
!			    deltaU = CurrentVelocities(j,CellOwner)-CurrentVelocities(j,CellNeighbour)
!			end if
!			do i=1,3
!				faceGradUij = meshSf(i,iFace) * deltaU
!				do k = 1,4
!					gradU(i,j,meshFaces(k,iFace)) = gradU(i,j,meshFaces(k,iFace)) + faceGradUij ! m2 * m/s
!				end do
!			end do
!		end do
!    end do
!
!	do iPoint = 1,nPoints
!	    if(PointViscousStencil(iPoint).eqv..false.) cycle
!
!		do j=1,3
!			do i=1,3
!				gradU(i,j,iPoint)=gradU(i,j,iPoint) / meshVPoints(iPoint) !  ! 1 / s
!			end do
!		end do
!    end do
!
!    do iFace = 1,nInternalFaces
!        CellOwner = meshFaceOwner(iFace)
!        CellNeighbour = meshFaceNeighbour(iFace)
!        if(nNextTimeCell1(CellOwner).eqv..false..or.nNextTimeCell1(CellNeighbour).eqv..false.)then
!
!            do j = 1,3
!                do i = 1,3
!                    faceGradU(i,j) = 0d0
!                    do k = 1,4
!                        faceGradU(i,j) = faceGradU(i,j) + gradU(i,j,meshFaces(k,iFace))
!                    end do
!                    faceGradU(i,j) = 0.25d0 * faceGradU(i,j)
!                end do
!            end do
!
!            if(nNextTimeCell1(CellOwner).eqv..false.)then
!                Rho = FaceVarsNew(1,1,iFace)
!                Temp = TemperaturePRho(FaceVarsNew(5,1,iFace),Rho)
!            else
!                Rho = FaceVarsNew(1,2,iFace)
!                Temp = TemperaturePRho(FaceVarsNew(5,2,iFace),Rho)
!            end if
!
!            divVel = faceGradU(1,1) + faceGradU(2,2) + faceGradU(3,3)
!
!            DynVisc = DynamicViscosity(Temp,Rho)
!		    BulkViscTotal = TotalBulkViscosity(DynVisc)
!
!            stress = 0d0
!            
!            do j = 1,3
!                do i = 1,3
!                    stress(i,j) = stress(i,j) + DynVisc * (faceGradU(i,j) + faceGradU(j,i))
!                    if(i.eq.j)then
!                        stress(i,j) = stress(i,j) + 2d0 * BulkViscTotal * divVel
!                    end if
!                end do
!            end do
!
!            do i = 1,3
!                ViscFlux(i) = 0d0
!                do j = 1,3
!                    ViscFlux(i) = ViscFlux(i) + stress(j,i) * meshSf(j,iFace)
!                end do
!            end do
!            
!            if(nNextTimeCell1(CellOwner).eqv..false.)then
!                do i = 1,3
!                    ViscousFluxes(i,1,iFace) = ViscFlux(i)
!                end do
!            end if
!
!            if(nNextTimeCell1(CellNeighbour).eqv..false.)then
!                do i = 1,3
!                    ViscousFluxes(i,2,iFace) = ViscFlux(i)
!                end do
!            end if
!
!        end if
!    end do
!    
!    ! external faces
!    
!    do iFace = nInternalFaces+1,nFaces
!        CellOwner = meshFaceOwner(iFace)
!        if(nNextTimeCell1(CellOwner).eqv..false.)then
!
!            do j = 1,3
!                do i = 1,3
!                    faceGradU(i,j) = 0d0
!                    do k = 1,4
!                        faceGradU(i,j) = faceGradU(i,j) + gradU(i,j,meshFaces(k,iFace))
!                    end do
!                    faceGradU(i,j) = 0.25d0 * faceGradU(i,j)
!                end do
!            end do
!
!            Rho = FaceVarsNew(1,1,iFace)
!            Temp = TemperaturePRho(FaceVarsNew(5,1,iFace),Rho)
!
!            divVel = faceGradU(1,1) + faceGradU(2,2) + faceGradU(3,3)
!
!            DynVisc = DynamicViscosity(Temp,Rho)
!		    BulkViscTotal = TotalBulkViscosity(DynVisc)
!
!            stress = 0d0
!            
!            do j = 1,3
!                do i = 1,3
!                    stress(i,j) = stress(i,j) + DynVisc * (faceGradU(i,j) + faceGradU(j,i))
!                    if(i.eq.j)then
!                        stress(i,j) = stress(i,j) + 2d0 * BulkViscTotal * divVel
!                    end if
!                end do
!            end do
!
!            do i = 1,3
!                ViscFlux(i) = 0d0
!                do j = 1,3
!                    ViscFlux(i) = ViscFlux(i) + stress(j,i) * meshSf(j,iFace)
!                end do
!            end do
!
!            do i = 1,3
!                ViscousFluxes(i,1,iFace) = -ViscFlux(i)
!            end do
!        end if
!    end do
!
!end subroutine CalcViscousFluxes
!
!
!subroutine DefineCurrentGlobalViscousStencil()
!use Data_Mesh, only: nCells, nFaces, nPoints, &
!                     nCell2NeigCells, Cell2NeigCells, &
!                     nFace2NeigCells, Face2NeigCells, &
!                     nPoint2NeigCells, Point2NeigCells, &
!                     CellViscousStencil, FaceViscousStencil, PointViscousStencil
!use DataTime, only: nNextTimeCell1
!implicit none
!
!    integer :: iCell,jCell,iFace,iPoint
!    integer :: i
!
!    ! compute cell global viscous stencil
!    do iCell = 1,nCells
!        CellViscousStencil(iCell) = .false.
!    end do
!
!    do iCell = 1,nCells
!        do i = 1,nCell2NeigCells(iCell)
!            jCell = Cell2NeigCells(i,iCell)
!            if(nNextTimeCell1(jCell)) cycle
!            CellViscousStencil(iCell) = .true.
!            exit
!        end do
!    end do
!
!    ! compute face global viscous stencil
!    do iFace = 1,nFaces
!        FaceViscousStencil(iFace) = .false.
!    end do
!
!    do iFace = 1,nFaces
!        do i = 1,nFace2NeigCells(iFace)
!            jCell = Face2NeigCells(i,iFace)
!            if(nNextTimeCell1(jCell))cycle
!            FaceViscousStencil(iFace) = .true.
!            exit
!        end do
!    end do
!    
!    
!    ! compute point global viscous stencil
!    do iPoint = 1,nPoints
!        PointViscousStencil(iPoint) = .false.
!    end do
!
!    do iPoint = 1,nPoints
!        do i = 1,nPoint2NeigCells(iPoint)
!            jCell = Point2NeigCells(i,iPoint)
!            if(nNextTimeCell1(jCell))cycle
!            PointViscousStencil(iPoint) = .true.
!            exit
!        end do
!    end do
!
!end subroutine DefineCurrentGlobalViscousStencil
!
!subroutine CalcExtrapolatedVelocities()
!use Data_mesh, only: nCells, CellViscousStencil
!use DataTime, only: nNextTimeCell, nNextTimeLocal, dnTimeCell
!use DataCells, only: CurrentVelocities
!use Data_Cfoam, only: Unew, U
!implicit none
!
!    integer :: iCell, i
!    real(8) :: alfa
!
!    do iCell = 1,nCells
!        if(CellViscousStencil(iCell).eqv..false.) cycle
!
!        do i = 1,3
!!            alfa = 2.0d0 * (nNextTimeLocal + dnTimeCell(iCell) - nNextTimeCell(iCell)) * dnTimeCell(iCell)**(-1)
!! linear extrapolation
!!            alfa = (nNextTimeLocal + dnTimeCell(iCell) - nNextTimeCell(iCell)) * dnTimeCell(iCell)**(-1)
!!            CurrentVelocities(i,iCell) = U(i,iCell) + (Unew(i,iCell)-U(i,iCell)) * alfa
!
!! Adam's method
!! U_n+1 = U_n + tau_n+1 * F_n = U_n + tau_n+1 / tau_n+1 * (U_n - U_n-1)
!            alfa = (nNextTimeLocal - nNextTimeCell(iCell)) * dnTimeCell(iCell)**(-1)
!            CurrentVelocities(i,iCell) = Unew(i,iCell) + (Unew(i,iCell)-U(i,iCell)) * alfa
!        end do
!    end do
!
!end subroutine CalcExtrapolatedVelocities
